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Abstract 

In this work we consider how surface-adherent bacterial biofilm communities respond in flowing systems. 
We simulate the fluid-structure interaction and separation process using the immersed boundary method. 
In these simulations we model and simulate different density and viscosity values of the biofllm than that of 
the surrounding fluid. The simulation also includes breakable springs connecting the bacteria in the biofllm. 
This allows the inclusion of erosion and detachment into the simulation. We use the incompressible Navier- 
Stokes (N-S) equations to describe the motion of the flowing fluid. We discretize the fluid equations using 
flnite differences and use a geometric multigrid method to solve the resulting equations at each time step. 
The use of multigrid is necessary because of the dramatically different densities and viscosities between the 
biofllm and the surrounding fluid. We investigate and simulate the model in both two and three dimensions. 

Our method differs from previous attempts of using IBM for modeling biofllm/flow interactions in the 
following ways: the density and viscosity of the biofllm can differ from the surrounding fluid, and the 
Lagrangian node locations correspond to experimentally measured bacterial cell locations from 3D images 
taken of Staphylococcus epidermidis in a biofllm. 

Keywords: Navier-Stokes equation, biofllm, immersed boundary method, computational fluid dynamics, 
multigrid, viscoelastic fluid 



1. Introduction 

In this paper we investigate the response and fragmentation of a biofilm attached to the interior of a 
tube and subjected to a flowing fluid. Specifically, we study the mechanisms of biofilm fiuid response and 
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detachment in terms of varying biofilm density, elasticity, and viscosity. In the simulations, we model biofilms 
attached to the walls of 3-dimensional square tubes using an extension of the immersed boundary method 
(IBM) (originally developed by Peskin Our approach differs from the traditional IBM in several ways. 

We use experimentally measured biofilm bacterial cell locations as initial positions for our Lagrangian nodes 
whereas traditional IBM refines the Lagrangian mesh along with the Eulerian mesh. As a result we also 
have to adapt the Dirac delta approximation to scale with the radius of the bacteria rather than with the 
mesh width. 

In this introduction, we first provide a brief background on the biology and biomechanics of bacterial 
biofilms in In Section L2 we discuss some alternative mathematical models that have been used to 



model biofilms (along with advantages and disadvantages). In Section 1.3 we introduce the immersed 



boundary method, and in |1.4| we discuss the significance of including variable viscosity. 

1.1. Bacterial Biofilms 

Biofilms are a phenotype of bacteria that are found in health, industrial and natural settings. In the 
medical field, biofilms occur on devices such as contact lenses, catheters, and mechanical heart valves. In 
industrial settings, they occur in and on water pipes, storage tanks, ship hulls, filters, food preparation 
facilities, etc. In natural settings, they can be found as slime on rocks in bodies of water or as dental plaque 
on teeth. 

Physically, biofilms are immobile and consist of a community of bacterial cells embedded in a dense 
surface-adherent extracellular matrix (ECM) of polysaccharides. Biofilms are mechanically strong structures 
that tend to deform and fragment rather than completely dislodge when subjected to fiows. Figure [T] on 
page [3] contains an electron micrograph of a biofilm of Klebsiella pneumoniae, clearly showing the ECM 
interconnecting the bacterial cells. The physical properties of the ECM are central to the growth, attachment, 
and detachment of biofilms. The focus of this work is on the biomechanical response of the ECM to fiuid 
fiow resulting in deformation and separation. 

An important feature of biofilms is that they are known to behave like viscoelastic fluids |17| . In other 
words, they exhibit both viscous and elastic responses upon deformation. Describing the exact viscoelastic 
behavior has been the subject of much experimental, theoretical, and computational research |17l 1551 |3J [551 
Uni Uni EU EH EH EZI • For example, Klapper et al. and Pavlosky et al. use a linear Jeffrey's constitutive law 
[171 while Lau et al. use a Voigt standard linear solid model for viscoelastic materials [20 . Our model 
includes elasticity of the biofilm by using simple linear springs (as done in [2 ) to connect the bacterial cells, 
and includes the viscosity of the biofilm with a modification of the constitutive equations for stress. 

1.2. Mathematical Models of Biofilms 

Much research into the mathematical modeling of biofilm growth and fluid/structure interactions has 
been conducted in the last three decades |33l[Mll5^H51imn51ll9| . Below, we summarize several modeling 
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Figure 1; Scanning electron microscopy of sessile K. pneumoniae LM21 performed in mature biofilm formed on Thermanox 
slides in microfermentor system after 48 liours of development at 20000 times magnification. This image is from Balestrino et 
al. [5] (used with permission). 

and simulation strategies. This is not intended to be exhaustive and we direct the interested reader to 
[441 116] for more in-depth reviews. 

The first attempts at mathematical modeling of biofilms were conducted in the early 1980s |36l 1551 [T5] . 
Picioreanu and others f l^ EU HH] advocated for an individual based (IB) approach, which models the 
behavior of each bacteria, encompassing ideas such as cell division, cell motility, metabolism, and death to 
simulate the growth and formation of colonies. Hybrid discrete-continuum models were the first methods to 
couple the flow with the the biofilm computationally in 2D and 3D simulations. Picioreanu, van Loodsdrecht, 
and Heijnen developed and used these hybrid discrete-continuum models to incorporate the flow over the 
irregular biofilm's surfaces, convective and diffusive mass transfer of substrate, bacterial growth, and biomass 
spreading [551 [Ml [52]. 

The most sophisticated (purely) continuum models developed are the phase-field models, which use a one- 
fiuid/two-component formulation in which the ECM and the bacteria are modeled as one fiuid component, 
while the collective ensemble of nutrient substrates and the surrounding fiuid are the other [461 147] . Two- 
dimensional simulations of both biofilm growth and biofilm-fiow interaction are presented in [47] , in which 
shear induced deformation and detachment are illustrated. 

We note that our model differs from both these approaches in the way that we model the biofilm. We treat 
the bacteria in the biofilm as discrete points, where the nodal locations in our simulations correspond to the 
locations of the bacterial cells within the biofilm. This contrasts from the continuum phase-field models that 
only include averaged biomechanical properties of the biofilm. With our mathematical formulation, just as 
in the individual based models, we can obtain the cumulative local stresses as well as attribute different local 
properties to the biofilm. Our model can be thought of as an extension of the individual based models, where 
we accurately account for the interactions with the fiuid as well as include the possibility of fragmentation. 
We also assume that on the time scale of our simulations there is no biofilm growth; so we ignore such 



factors as nutrient concentrations and growth rates. 

1.3. Immersed Boundary Method 

The overall goal is to simulate the response of a biofilm attached to the walls of both 2- and 3-dimensional 
square tubes. We do this using an extension of the immersed boundary method. In this section, we introduce 
the immersed boundary method and provide the framework necessary for us to later extend the IBM for 
our use. 

The immersed boundary method (IBM) was originally developed by Peskin to study blood flow in the 
heart [2^. The IBM has been used previously to model and simulate biofllm/fluid interactions by Dillon, 
Fauci, et al. in [12 , and by Alpkvist and Klapper in [2J. The authors successfully coupled the fluid to 
the biofllm; however, they make the assumption that the biofllm has the same density and viscosity as the 
surrounding fluid. This choice substantially simplifles the task of solving the N-S equations but does not 
account for the fact that biofllms typically have 500 x larger viscosity and 12% larger density than water 
[171 137j . They also use a random distribution of points within a biofilm-shaped shell to represent the biofilm, 
which does not account for the true spatial distribution of bacterial cells within a biofilm. 

The IBM has been used more recently in the modeling of immersed elastic structures in viscous flows 
in [HI W^ . in which the authors use constitutive viscoelastic models including Maxwell, Voigt, and 

Jeffrey's models to incorporate forces in the immersed structures into the IBM. Similarly, our ultimate goal 
is to establish an appropriate constitutive model for the forces in the biofilm with the help of experimental 
collaborators and to include this in our IBM formulation. 

1.4. Variable Viscosity 

It is agreed upon in the biofilm research community that biofilms behave like a viscoelastic fiuid. There 
has been a few efforts to match the behavior of biofilms with conventional mechanical viscoelastic models 
|17 [ l3| [28 | l20|. However, these efforts have not produced a consensus on the model to use for viscoelasticity 
in biofilms. This is, in part, due to the fact that the viscoelastic properties in biofilms is highly variable 
with different growth conditions |10| and even in the same growth conditions [T]. 

The incorporation of spatially variable viscosity in the immersed boundary method is an area that has yet 
to be well developed. Luo et al. couple the immersed viscoelastic structure to the fluid flow in an immersed 
boundary type formulation, but they solve the fluid equations separately from the equations governing the 
motion of the immersed viscoelastic solid and then couple the solutions at their physical interface [23] . This 
formulation will not work in our case because we couple the biofilm to the fiuid in the entire domain so that 
it will behave as a viscoelastic fiuid. Another approach to including viscosity into the immersed boundary 
method is by replacing the simple elastic springs with viscoelastic links, which will change the value of the 



external force, f, in the Navier-Stokes equations (Equation (26 1 below). This type of strategy was used first 
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by Bottino in [6 to model general viscoelastic connections in actin cytoskeleton of ameboid cells and also 
by Dillon and Zhuo in |49| to model sperm motility. 

There are two natural ways to add viscosity to the biofilm. The first is through the use of local damping 
forces in addition to the springs that we use for the elastic component. In this way, we can define the forces 
between any two connected bacterial cells using a typical mechanical model for a viscoelastic material. The 
second way is to treat the entire domain as a continuous Newtonian viscous fluid with a spatially varying 
viscous coefficient. The use of dashpot damping with our current mathematical formulation is quicker to 
implement but has serious stability restrictions in the simulation, and thus we omit this method herej^ 

The approach in this work is to treat the fluid in the entire domain as a Newtonian viscous fluid with a 
spatially varying viscous coefficient. The core idea involves replacing the viscous term in the Navier-Stokes 
equation with one that can apply to nonuniform viscosity in the fluid. We note that to the best of our 
knowledge this approach has not yet been attempted in the single fluid immersed boundary method. There 
have been attempts using two materials (fluid-fluid or fluid-solid), coupling them at their interface, in which 
the stress is adapted in the viscoelastic fluid or solid to account for a different viscosity |47l However, 
in our approach, we couple the biofilm to the fiuid within the entire biofilm region, not just at the interface. 
Thus, we must adapt the forces in the Navier-Stokes equations to account for the viscous and elastic stresses 
on the surface and within the biofilm. 

We now describe the organization of this paper. In Section [2] we provide our mathematical formulation, 
which is a variation of the immersed boundary method literature. In Section |3] we describe our numerical 
method, based on a multigrid approach. In Section |4] we provide numerical validation of our method. In 
Section [5] we provide simulation results in both two and three dimensions, running our simulations for 
a variety of experimentally obtained biofilms with varying parameters such as spring constants, densities, 
viscosities. Finally, we provide conclusions in Section [7] and a discussion of possibilities for future work in 
Section [i 

2. Mathematical Formulation 

In this section, we provide the mathematical formulation for our simulations. We use an Eulerian mesh 
to describe the system as a whole and solve the dimensionless N-S equations at each time step on this 
mesh. The Lagrangian nodes are used only to compute information about the biofilm (location, velocity, 
local density, force) and then transfer the information back onto the Eulerian mesh using the Dirac delta 
function For convenience, we provide a list in Appendix A of the variables and parameters used in this 
work. 

^We will pursue this approach in future work when change our numerical scheme to a se mi- implicit or implicit method. 
^The Dirac Delta function is approximated in the actual implementation (see Equation Jl4|). 
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We now introduce the mathematical equations used in our model. The dependent Eulerian variables 
are velocity u(x, t), pressure p(x, t), density /9(x, t), and Eulerian force density f(x, t), where x is the 
independent Eulerian variable and t is time. The dependent Lagrangian variables are position of the nodes 
X(q, t), velocity of the nodes U(q, t), and the Lagrangian force density F(q, t), where q = {q,r,s) is the 
independent Lagrangian variable. The equations of motion for the biofilm-fluid interaction are 

p(x, (^^+u-Vuj = -Vp+V- (/i(x, ^)(vu+(Vu)^)) +f(x, t) (1) 
V • u = 0, (2) 
— (q,t) = U(X(q,t),i), (3) 

f(x, = / F(q, t)5(x~X(q, t))dq, (4) 

p{x,t) = po+ ph(5(x- X(q, t))dq, (5) 

U(X(q, 0,0 = /ulx, 0<5(x-X(q, 0)dx, (6) 
Jn 

where /i is the dynamic viscosity, po is the mass density of the fluid, pi, is the additional mass density of the 
biofilm from that of the surrounding fluid, is the flow domain, fif, C fi is the space occupied by only the 
biofilm, and 5{x) is the Dirac delta function. Equations ([T]) and ^ are the incompressible Navier-Stokes 
(N-S) equations with spatially varying viscosity and a forcing term that represents the forces applied by the 
biofilm on the fluid. Equation (jsj is the equation of motion of the biofilm, where U(q, t) is the velocity of 
the biofilm. The systems of PDE's given by ([l])-(|2| is coupled to ([3| by the integrals given in (|4|-([6]). 

To avoid numerical inaccuracies due to roundoff errors, we non-dimensionalize these equations using the 
non-dimensional variables defined as 

f* _t_ Y* 2i 11* — _iL n* P~P'^tv.ba 

Po ' fa' ^ Mo ' 

where po is the pressure at the upstream end of the tube, PLt^h^ pressure at the downstream end 

of the tube, T is the characteristic time scale, /o is the characteristic force density, and L is the character- 



istic length. We use the scaling parameters defined in Table [BTTO] on page|43) Dropping the stars from the 
dimensionless variables, equations ^ and (|4])-(|6]) remain the same as in the case with dimensions, while 
equations ([T|) and ([3| become 

crp(x, 0^ +P(x, t)u- Vu = -eVp + i^e^^V • f^(x, <) fvu+ (Vu)^)) + ^^f(x, i), (J) 
ut V V / / Po^o 

a — (q,t) = U(q,i), (8) 

where a — is the Strouhal number, e — ^° ^^i^&g jg the Euler number, and Re = p°^"° is the Reynolds 
number of the fiuid. 

The initial velocity profile is the exact solution to the incompressible Navier-Stokes equations in a square 
or circular tube with rigid walls and no-slip conditions at the walls. The velocity profile for a circular cylinder 
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can be found in many textbooks in fluid dynamics (such as [45j), and a series solution for the laminar flow 
velocity profile for a square tube was derived by Spiga and Morini in |40,. 



3. Numerical Method 

In this section, we describe the numerical formulation for our simulations. Our numerical task is to solve 
the system defined by Equations ([T])-([6]), and we now provide the details of our numerical approach. 

The incompressible fiow Navier-Stokes equations, ([lj-(|2|, are discretized on a fixed uniform Eulerian 
lattice, while the biofilm equations are discretized on a moving Lagrangian array of points that do not nec- 
essarily coincide with the fixed Eulerian mesh points of the fiuid computation. We represent the interaction 



equations (|4|-(|6| with a smoothed approximation 5 to the Dirac delta function (see 3.11. Our numerical 
approach was inspired by the solving technique used by Zhu and Peskin in [48J to simulate a fiapping filament 
in a soap film. 

The discretized equations corresponding to (|4|-([6| are given by 

r(x) = f]F"(s)5~(x-X"(s),..), (9) 

s=l 

p"(x) = po + J2PbSi^-^"i-^)^^)dl (10) 
U"+i(s) = ^u"+i(x)^(x-X"(s),/i)/i3, (11) 

X 

where the superscript n denotes numerical approximations at a particular time step n, rj is the total number 



of Lagrangian discretization points, the sum in ( 11 1 is over all the discrete points of the form x = (i/i, jTi, kh) 
with i, j, and k are integers, h is the Eulerian mesh width, and df, is the average volume element of the 
Lagrangian nodes (computed by dividing the total volume of the biofilm by the total number of Lagrangian 
nodes distributed within it). Following convention, we replace (q, r, s) from the mathematical formulation 
with only s, which we use as an indexed label with a unique number assigned to each Lagrangian point |48| . 
In (|9|, F(s) is now the total elastic force on the Lagrangian node associated with marker s, as opposed to 
an elastic force density. This is because we calculate the force explicitly depending on which other nodes it 
is connected to. 

3.1. Dirac Delta Approximation 
In [30], Peskin defines (5/i(x) as 

«.,x) = k-V(^)*(|)*(|) , (12) 



where 4'{r) is 



if Irl < 1, 



|(3-2|r| + Vl+4|r|-4r2); 
(r) = (/.i = <j i (5-2|r|- v/-7+12|r|-4r2) . if 1 < |r| < 2 , . (13) 

0; if|r|>2. 
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We replace this Sh, that is used in standard IBM implementations, with one that scales with uj instead 
of h as 



(14) 



We deviate from the standard scaling of the Dirac Delta approximation for two reasons. The first is that 
we wish to give a presence to the bacterial cells that is representative of the true volume of the cells. Thus, 
in the simulations, we make uj in ([9]) and (10) equal to the radius of a bacterial cell that we are modelling. 
Equation ^ then spreads the force over a volume that is slightly larger than the cell, ensuring that the 
entire space occupied by the cell in the fluid is influenced by the force. The second reason we use this scaling 
is because, during the mesh refinement analysis described in |4.3.3| we discovered that the implementation 
with the scaling by h restricts us to less than first-order convergence of the velocity, u. Using a scaling that 
is independent of the mesh-width fixes this issue and leads to greater than first order convergence. 

However, this modification does have the negative consequence of losing two desirable conditions that 



were previously satisfied by the Dirac delta approximation, 6^- Specifically, with as defined in (13 1, the 
unity condition, 

^ 5,,(x-X)/i3 = l, VX, (15) 



and the first-moment condition, 



^ (x - X)(5,,(x - X)h^ = 0, VX , 



(16) 



are both satisfied. However, using S in place of Sh, these conditions fail to hold true for all X when co ^ h. 
In practice, this is not a major concern as many IBM formulations use a Dirac delta approximation that 



satisfies the unity condition, (15 1, but does not satisfy the first-moment condition, (16 1. For example, in 



], Peskin and Zhu replace 4>{r) in with 



i(l + cos(f)); if|r|<2, 
0; if|r|>2. 



(17) 



We do have to choose whether we want to use ^^(r) as defined by ( 13 1 or by ( 17 1. For either choice of 0, it 
is true that both 

' S{x - X)dx = 1 



lim 



and 



lim 

/i-i-O 



(x ~ X)(5(x - X)<ix = . 

Note that, in the limit as /i — 0, we see greater than 0{h^) convergence to (15 1 and (16 1 (see Figure [2] on 
page 10 1, which is consistent with the theoretical convergence rate for a Riemann sum. Therefore, we choose 



to use the for which the summations in (15 1 and (16 1 are closest to 1 and 0, respectively, for the values of 



Lo and h used in our simulations so that we can have the most accurate discrete approximation of the Dirac 
delta function. 

We now define two error metrics to determine how well 5 (using either 0i or (1)2) satisfies the unity and 
first-moment conditions. Analogous metrics and comparisons could be conducted in higher dimensions, but 
for simplicity we provide a one-dimensional comparison. Using cj = and the one-dimensional version of 



(14 1, we define 

" (5(x-X, uj)h\-l 



(18) 



— X) (5(x — X, uj)h 



(19) 



e«mty(w, h) = max 

and 

tmom{LO,h) = max 

We only have to find the maximum over X G [0, ft.] since the summations are periodic with period h, because 
X = and X = ft both correspond to a Lagrangian point being at the same location as an Eulerian point. 
We find eu7iity{^/^0Q, h) and e„io„i(Yioo, ft) for values of ft e {jf^i ifjo) ^^id using both (pi and 02. These 
values are compared in Figure [2] on page |10| and show that using 02 provides a better approximation of the 
Dirac delta function in terms of matching the values of these summations for most values of ft. Notice too, 
that the values of Cunity and emom are exactly zero for certain relationships between ft and lo. For example, 
if ft = ^ with z G N, then tunity and tmom are exactly zero when using 01. Thus if our problem allows 
for h — J, then we would use 0i. However, with the choices we have made for uj and ft in our numerical 
scheme, it is a better choice to use 02 in our simulations. 

In the transfer of information from the Lagrangian grid to the Eulerian grid, we scale by a; (see (|9| 



and (10)). However, in the transfer of velocity from the Eulerian to Lagrangian grid (Eq. ( |11[ )), we scale 
with ft instead of oj in order to capture the velocity only at the center of mass of the bacterial cells. 

3.2. Elastic Forces 

In [2], Alpkvist and Klapper use Hooke's Law to describe the elastic force between the connected La- 
grangian nodes. We also use this method as our first attempt to model the interconnecting force in the 
biofilm. Thus, the elastic force on each Lagrangian point using Hooke's Law is 

r"(s) = ^/,,fc^T,,fc, (20) 

k=l 

where T is the tension between nodes s and fc, / is the connectivity matrix defined as 

1 bacteria s connected to bacteria k 
otherwise , 

and ds^k is the vector pointing from Lagrangian node s to A: with magnitude dg^k- The tension from the 
spring connecting node s and k is formulated as 

Ts,k = Ks^k{ds,k — ^s^k)} 

9 



I0g2(/J) 

Figu re 2: We show ^unity (i/lOO, h) and emomiy^OO, h) for h e (j^, j^jo)- '/'I t^e </> given in |l3| l and (/)2 is the 4> given in 
\n\ . We show log2 in the x and j/ axes so that the convergence rate appears as the slope of the line segments. 

where r^- fe is the rest length of the spring connecting nodes s and fc, and K^^k is its Hookean spring coefficient. 
We choose to define each spring coefficient as 

= — , (21) 

rs,k 

where Fmax is the force required to brealc the spring. We define the spring coefficients in this way to ensure 
that all of the springs, regardless of initial length, break with a force of Fmax when they are stretched to 
a length of 2r<;_fc. In our simulations, we vary F„iax to attain specific results (such as detachment; see 



5.1 



more details). As is done in [2J, we model the failure of the ECM by breaking the connections between the 
Lagrangian nodes as the springs used to connect them exceed twice their resting length. We note that this 
condition, however, is not based on experimental evidence, and in future work we will adapt this breaking 
criteria according to experimental results. In Section [8] we discuss future adaptations to the breaking criteria 
in terms of the yield stress of polymers. 

We conclude this section with a ID illustration to show how 3 linearly connected cells transfer their 
elastic forces to the Eulerian grid via Equation (|9| , and how this is effected by varying h (see Figure |3] 
on page 111. In this example, the three nodes are close enough to each other that their forces add in the 
overlapping regions. Sub-figures (b)-(e) illustrate how, with finer discretizations (smaller /i), the elastic 
forces on the bacteria are more accurately represented in the Eulerian grid. 
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Figure 3: Examples of ID f{x) with three Lagrangian nodes (a), one at X{s = 1) = ^/32 on an Eulerian node and the others at 
X{s = 2) = ^ + and X(s = 3) = ^ + ^io (this position is chosen so that the second and third nodes are close enough to 
the first one to show the interaction of the forces), with u) = These positions were chosen so that the cells would be close 
enough to each other to have an overlapping region on the Eulerian grid after the transfer of the forces from the Lagrangian 
grid. In this demonstration, F{s = 1) = 5, F{s = 3) = -8, and F{s = 2) = -(F(l) + F(3)) = 3. These plots illustrate the 



effect of using different spatial steps; (b) h = (c) h = -gj, (d) h ■■ 



{e)h: 



1 

512 • 
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(a) 
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(c) 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



(d) 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 



Figure 4; Examples of ID fi{x) with two Lagrangian nodes, one a,t x = 5/32 on an Eulerian node and the other at x = + j^, 
with LU = These Lagrangian nodes (cells) are close enough that there is a region of interaction. These plots illustrate the 
effect of using different spatial steps: (a) h = (b) h = ^ (c) /i = j^g (d) h = 



3.3. Variable Viscosity 

It is known that ECM density decreases with distance from an individual cell. To account for this, the 
exact form of /i(x) used in our simulations is 



/i(x) 



max 

l<s<r), sGP 



fJ.out)S{x - X(s), Uj) + Hout 



(22) 



where ^max is the viscosity at a bacterial node, ^out is the viscosity of the surrounding fluid, D is the spatial 
dimension, and a; is a parameter we can use to stretch the influence of the additional viscosity. We made 
this choice for /i(x) because we wanted a viscosity that would decrease at the same rate as the elastic force 
with the distance from the bacterial cell. See Figure |4] on page [12] for a ID example of the effect of lo and h 
on the viscosity distribution, /x(x), from two interacting cells. In the future, we can change this function to 
suit the specific viscous properties of the particular biofilm. 
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3.4- Solution Strategy 

We employ a projection method to solve the incompressible Navier-Stokes equations numerically, 
building on the method used by Zhu and Peskin in [JS]. This method introduces a velocity field (at an 
intermediate time). u(x, t), which is the solution to the difference equation 

(" "'^IP' + ^" ■ + ■ (""^-))") ^ • + Dl.n--^')) + ^f," , (23) 

where k — 1, 2, 3 and the subscripts k denote the k^^ component of that vector. The finite difference 



operators in ( 23 1 are originally defined by 



(x + hei) + (f)(x — he-i) — 2(/)(x) 



DIM = 



(x + hei) — (/)(x — hei) 



2h 

— \^h,n ^h,2T ^h,3)y 



where is the unit vector in the i*'' direction. Additionally, D° • (aD^i/)) is defined for scalar functions, 
a(x) and 4>{x.), using the midpoint values of a as 



D".(aD»=^ 

1=1 

and D" • ^al?^ ^.u^ is defined as 
D°-(ai?^u) = 



^ a(x+|e,: 



(x+/iei)-0(x) 
h 



- a (x - |ei) 



0(x) — (/i (x— /iGi ) 



(24) 



El — 



2/i 



(25) 



2/i y 

Note that in the case of constant viscosity, using the incompressibility constraint, ([T]) is reduced to the 
standard N-S equation used in IBM, 



i) ( ^ + u- Vu 



-Vp + ^Au + f(x, t). 



Therefore, in the case of constant viscosity, equation (23 1 is replaced with 



P o- 



~Tt+l n 



At 



— + DOufc + D" • (uufc))" j = ^^Lh{n 



.~.n+l\ 

k ) 



Lfo en 
^li- ■ 



(26) 



(27) 



We solve ( 27 1 as opposed to ( 23 1 in the constant viscosity case, since each component of u, , is independent 



of each other in (27 1, and can therefore be solved for separately and simultaneously. 
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To complete the discretized incompressible Navier-Stokes system, we have the following two equations: 

= -eD V+i , (28) 



up 



At 

DO-u"+i - 



0. 



(29) 



We point out here that summing equations ( 23 1 and ( 28 1 leads to the discretized version of (261, with the 
exception that the evaluation of the viscous term is at the intermediate value of the velocity. We solve for 



pressure by applying D*^ to both sides of ( 28 1 and using ( 29 1 to obtain 

DO • ( — DV 



P" 



a D" • u"+i 



(30) 



To solve equations (23 1 and (30 1, we use Gauss-Seidel as a smoother in a multigrid solver. At each time 
step, we solve (27 1 for u"+^, substitute it into (30 1, solve for and finally solve for using (28 1. 



Then the velocity is transferred from the Eulerian points to the Lagrangian points using (11 1. With U"+^ 
computed, the new Lagrangian node locations are computed using Euler's method as 

X"+i(.s) = — U"+i(.s) +X"(s). 
a 

The forces between the Lagrangian points are then recalculated and transferred to the Eulerian points using 
([9|. Finally, the values of p" and /x" are evaluated using the new Lagrangian locations. 

3.5. Multigrid 

In this section, we discuss the elements of multigrid that we use in our solution strategy. For more details 
on multigrid, see [7]. 

In our solver, we use the conventional Gauss-Seidel iterative method with red-black ordering. In the 
multigrid scheme, we use full-weighting restriction to go from fine to coarse grids, and we use linear interpo- 
lation to go from coarse back to fine grids. The finest grid is the grid with step size h and the grids become 
coarser by increasing the step size by a factor of 2. This halves the number of nodes in each dimension, 
allowing for significantly faster computations on the coarser grids. The number of levels in the multigrid 
solver depends on both the dimensions of the computational domain as well as h. In our simulations, we 
iterate using multigrid V-cycles until we reach a sufficiently low value for the norm of the residual, 

at each time step. Here, A'^v''^ — is the linear discretization of a PDE, and the residual is r'' — f''' — A^v'^, 
where v is an approximation to v. The residual provides a bound on the true error in the solution of the 
linear system since we have this relationship between the error and the relative residual error: 



< cond{A'' 



in ' 



(31) 



where e 



_ jh _ j[h~h^ g^^^ cond{A'^) is the condition number of A'^ |T. In Section [J] we give 



approximations for the condition numbers for our matrices. 
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(a) ' (b) 



50. 

45^ 10- 




(c) . . (d) 

Figure 5: Shows examples of the computational domains with a sample biofilm. Axes units are microns, (a) is 2 dimensional 
and (c) is 3 dimensional, (b) and (d) show enlarged images of the biofilms. The points shown in these plots show the initial 
position of all of the bacterial cells in the biofilm. 



3.6. Boundary Conditions 

The computational domain used in our example simulations is a section of a tube with the biofilm 
centered in the direction along the axis of the tube (see Figure [s] on page 15 1. In the 2D case, flow is along 
the X-axis and, in 3D, it is along the z-axis. The boundary conditions we used in these simulations were 
derived from exact solutions for the velocity and pressure in the case of laminar flow. We now provide the 
boundary conditions in both the 2D and 3D cases. 

3.6.1. 2D Boundaries 

The no-slip boundary condition exists at the walls of tube and requires that the velocity be zero there, 
so we use that as the boundary condition at the walls. The velocity at the upstream boundary (x = 0) is 
held at the laminar flow velocity (shown in Figure [6] on page |17[ a)) given by 

ui{y) = ^{y^-2ay), (32) 

where a is the radius of the tube, y is the displacement from the bottom edge of the tube, k is the linear 
rate at which the pressure decreases through the tube, and ui is the x-component of the velocity (i.e., 
u = (ui, U2)). At the downstream boundary, a Neumann condition is applied to the velocity by enforcing 
that 

d \ 

where Xdown represents the x value at the downstream boundary. 
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The boundary conditions for pressure come from the laminar flow equation for pressure given by 



p{x) = Kx+p{0). (33) 
In our simulations, we hold the pressure at the upstream boundary at p(0) and at p{xj^own) the down- 



stream. At the top boundary, we hold the pressure at the values given by (33 1 and, at the bottom boundary 
(the boundary on which the biofilm is attached) , we use a Neumann boundary 

d \ 

9y Jy=0 

3.6.2. 3D Boundaries 

In the 3D simulations, we orient the square tube along the z-axis (see Figure [s] on page [Ts^c)). The 
no-slip boundary condition exists at the walls of tube and requires that the velocity be zero there so we 
use that as the boundary condition at the walls. Derived by Spiga in j40) . the velocity at the upstream 
boundary is held at the laminar flow velocity (shown in Figure [6] on page |17[ b)) given by 

IQko^ v-^ sin (mr^/a) sin (mirv/a) 
Mx,y) = ^ 2^ ^ 1 2^ ' (34) 

n,m>0,odd ^ ' 

where a is the width of the tube and is the z-component of the velocity (i.e., u — (ui, U2, us))- At the 
downstream boundary, a Neumann condition is applied to the velocity by enforcing that 

d 



u(x,y,z) = Vx, y. 

The boundary conditions for pressure come from the laminar flow equation for pressure given by 

p(z) = Kz+p(0), (35) 

where z = is the upstream boundary. In our simulations, we hold the pressure at the upstream boundary at 
p(0) and at p{zdown) at the downstream. At the top and side boundaries, we hold the pressure at the values 



given by Equation (35 1 and, at the bottom boundary (side with attached biofilm), we use the Neumann 
condition given by 



d 

p(x, y, z) ) =0 Vx, z. 

y=0 



dy 



4. Validation 

It is crucial to validate our numerical method using known results and mesh refinement convergence 
analysis. Thus, below we provide numerical evidence that both our 2D and 3D simulations and numerical 
methods are working as they should. 
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(a) (b) 

Figure 6: We illustrate upstream boundary conditions: (a) 2D laminar flow velocity profile, (b) 3D laminar flow velocity proflle. 
Units are in meters for (a) and (b). 

For the purposes of the 2D and 3D convergence analysis conducted in later sections, we require the 
following notation. We present the following notation in 3D (the 2D versions are analogous but with- 
out the z elements). Define the Eulerian grid function p-norm for an arbitrary 3D vector field, w(x) = 
(w;i(x),W2(x), W3(x)), by 

IIHIp= (EI^(^-%'^'»)I''^'') ' (36) 
where D is the spatial dimension, 1 < P < oo, and 

\w{xi,yj,Zk)\ = ^wi{xi,yj,Zky + W2{xi,yj, Zk)^ + W3{xi,yj, Zk^. 

Then 

ll^lloo = max\w{x„yj,Zk)\ . 

Additionally, on the Lagrangian grid define the Lagrangian grid Junction p-norm for a vector field, X — 
(Xi(s),X2(s),X3(s)), as 

ii^iip= (Ei(^i(*)'^2(s),x3(s))r<j , 

where 1 < p < oo and is the average volume element of the Lagrangian nodes. Then 

||X||^= max |(Xi(s),X2(s),X3(s))|. 

l<s<?7 

Note that both of these grid function norms are derived from using discretizations of the integrals used in 
a typical function p-norm (see Appendix A of [21J for more details). 

There are three parts to our simulation validation process: 1) we illustrate that in the absence of 
the biofilm our numerical simulation converges to the analytical solution; 2) we verify that our multigrid 
technique is correctly accelerating the convergence of our chosen relaxation scheme; and 3) we determine 
the convergence rate of the simulations with a biofilm using a mesh refinement convergence analysis. 

Before discussing the results of our validation process, we provide a brief description of the two primary 
sources of error present in our simulations, |4.1[ We also setup our simulations with a detailed description of 
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initial Lagrangian node positions in |4.2[ Then we provide 2D simulation validation in |4.3| and 3D validation 
inl44l 



4-.1. Discussion of Errors 

In our numerical scheme, we have two sources of error: 1) discretization error is introduced by discretizing 
the Navier-Stokes equations in space and time; and 2) algebraic error is introduced when we attempt to 
solve the resultant systems of linearized equations. 

As it is impossible to compute the true algebraic error, we use the norm of the residual to deduce an 



upper bound on the algebraic error using Equation (31 1. Recall that Equation (31 1 indicates that the 
relative algebraic error at each timestep is no larger than the condition number of the matrix times the 
relative residual norm (recall that jjj-^ is the relative residual norm). We do not construct these matrices 
during the actual simulations because we do not need them to solve the systems. However, we did construct 
them to find their condition numbers and found that the condition numbers for the matrices used in the 
mputations for u"+^ and p"+^ are O (h^^) (this is true for both 2D and 3D simulations). The simulations 



co: 



that resulted in the plots given in 5.2 and 5.3 were run using h = j^, and thus the matrix condition numbers 
were approximately 10^. 

Our goal in the simulations is to ensure that the algebraic error falls well below the discretization 
error at each time step, so the total error will be dominated by the discretization error. In theory, the 
discretization error is at best 0{h^) w C (i/iss) w C x 6 x 10"^ for C > 0, with our discretization. Using 
a stopping criteria of 10^^ for the relative residual at each timestep should suffice (i.e., from (31 1). We 
continue to the next time step only when the computed relative residual, 
lie'' II < cond{A) * IQ-^ w 10"^ 



< 10 ^, because this implies 



Another factor influencing the capability of our simulations is that after, extensive simulation, we discover 
that our linear solver is limited to converging to a relative residual norm of about 10^^^ (possibly from 
machine precision issues). With h — ^j^^j the condition number is O(IO^) and the discretization error is 
O(10~^), so the algebraic error is at best bounded by about 10^^^ x 10^ — 10^® (see Equation (31)), and we 
can no longer be certain that the algebraic error falls below the discretization error at each timestep. For 
this reason, we restrict h to be larger than in all of the simulations and convergence analysis. 

4-2. Simulation Setup 

In these convergence simulations, we constructed an experimentally motivated mushroom shaped biofllm 
(shown in Figure[7]on page |19[ a)). We carved this shape from a 1.6 fim slice cut from data points generated 
in the Younger and Solomon labs at the University of Michigan. These data points are 3D bacterial cell 
locations from 3D Leica SP2 confocal laser scanning microscopy images taken of Staphylococcus epidermidis 
RP62A (ATCC 35984) grown in a Stovall 3 channel flow cell for 24 hours at 37 C under a wall shear stress 
of 0.01 Pa. For further details of how the coordinates were computed, see Stewart et al. |41| . 
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(a) " (b) X ^ 

Figure 7: Mushroom-shaped biofilm at t = in the middle of the computational domain attached to the bottom (y = 0) of the 
tube, (a) This is the shape used in the 2D simulations, and (b) 3D mushroom shaped biofilm at t = 0. 

From this data, the average Lagrangian volume element, dp. is calculated to be approximately 4.036 finr^, 
and thus we use do = 1.59 /im in both the 2D and 3D simulations. We connect the initial distribution of cells 
with a distance based connection criteria. Our inspiration for the connection distance criteria came from 
the closeup images of biofilms such as the one shown in Figure [T] on page[3j We observed that each bacterial 
cell is connected to neighboring cells that are within about 2do- Thus, we varied the connection criteria in 
our algorithm between 1.5-2.5 x do in an effort to find one that resulted in a biofilm that was sufficiently 
connected but not overcrowded. This resulted in the choice of a connection criteria of dc — 2.8 /im. In other 
words, we placed spring connections between Lagrangian nodes at the beginning of the simulation with 
every node connected to every other node less than 2.8 /xm away. Admittedly, this value of d^ is arbitrary, 
and future work will include deriving a method to determine this connection criteria through image analysis 
of closeup images of biofilms similar to Figure [T] on page |3] The mushroom shaped biofilm has a height 
of about 8.5 /xm and width of about 8/xm (see Figure [t] on page 19 1. In the convergence simulations, the 



maximum spring force, F„iax, is set to 5 x 10^^ N. The fluid parameters for these convergence simulations 
are provided later in Table [9] on page [28] 

Note that (5 is a function of w, a scaling parameter we must choose that determines the volume/area of 
influence when the forces and density are transferred to the Eulerian grid. We must also point out here 
that a more accurate representation of the Dirac delta function occurs when lu > h. Thus, for the purpose 



of the convergence simulations, we use uj = 1.0 /im in the transfer equations, ([9| and (10 I. However, in our 
simulations, we use uj — 0.5 /im since the actual radius of Staphylococcus epidermidis is known to be about 
0.5 /im [43_. Using a characteristic length of L = 50 /im, we have the non-dimensionalized oj* = ^ = 
Dropping the star from the dimensionless variable, we use w = ^ in the convergence simulations and 
UJ = in the results simulations. We desire that u > h, so that the Lagrangian forces are spread at least 
two Eulerian mesh widths in every direction (as is done in the traditional IBM |30j). Using ui = in the 
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Figure 8: Exact error in biofilm-free 2D simulation with h = The solution is compared to the exact solution, |32| , using 

a maximum norm. 

convergence simulations allows h ^ , to obey these criteria. We again note that one of the reasons 



for using u! in the scaling of (14) as opposed to h is that better spatial convergence rates are achieved since 



the scaling is independent of h. 

4-.3. Two-Dimensional Validation 

In the absence of a biofilm, we expect that using a centered finite difference approximation for the second 



derivatives allows exact convergence to the second-order polynomial solution (Equation (32l). That is to 



say, we expect the numerical solution to converge to the analytical within machine precision. Reassuringly, 
we find that the biofilm-free simulations converge exactly to the steady state laminar fiowj^To illustrate, 
we started with an initial velocity profile that is one-half of that of the laminar fiow velocity profile. The 
error in the simulation converged (within machine precision) in less than 300 time iterations for all spatial 
resolutions (see Figure [s] on page 20 for example with h = 1/128 and dt — 0.0001). 

4- .3.1. Multigrid performance 

Next, we provide numerical evidence that the multigrid technique convergences optimally to the solutions 



of (27) and ( |30[ ). Define a work unit as the cost of performing one relaxation on the finest grid (see |7])- 
In Figure |9] on page 21 'a), we depict (for the pressure computation) the work units required to reach the 
minimum residual error as a function of allowed levels in the multigrid. This result shows that the number 
of work units required decreases significantly with each added multigrid level. This means that the multigrid 



*In the 3D simulations, we do not see exact convergence since the laminar solution is not a second-order polynomial. See 
|5.3| for details on the convergence properties of the 3D laminar flow case. 
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method correctly accelerates the convergence of our iterative method by doing computations on the coarser 
grids. For example, with just one allowed level of multigrid, the relaxation uses only the finest resolution 
grid and requires about 10^ work units, whereas with 6 multigrid levels we only require about 10^ work 
units to achieve the same error. Note that there is no reduction in the number of required work units with 
the addition of a 7*'^ level in the multigrid, so we use at most 6 levels in our 2D solvers. The data in this 
plot was obtained using our 2D simulation with a mushroom shaped biofilm similar to those shown in |5.2| 



O 
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(a) 



Number of Levels in the Multigrid 



(b) 



Number of Levels in the Multigrid 



Figure 9: Decrease in work units required to reach the minimum residual error as the number of multigrid levels is increased 
in the (a) 2D simulations, (b) and 3D simulations. 



4-3.2. Empirical Estimate of Convergence Rate in Time 

Similar to the development by Mori and Peskin in [25] , we define a measure of error by 

EMn^t)^ q^\T) ~ q^'^^T) , (37) 

p 

which is the error difference at time t — T in a, computed quantity, q, using a temporal refinement of a half 
timestep. Then, an empirical estimate for the convergence rate is calculated using 



EpiqiT);At) 



rMn^t) = log2 ( JZJ,[J, I . (38) 



We compute the approximate convergence rate in time using the E2 and Eao errors in the Eulerian variable, 
u, and in the Lagrangian variable, X. We simulate until t = T = 0.01 s using temporal step sizes that ranged 
from A< = 1/5000 to At = i/soooo, decreasing by a factor of 2 at each level. The Eulerian grid is discretized 
with a step size oi h = 1/256. 

The empirical convergence rates from our temporal refinement are provided in Table [T] on page |23[ The 
immersed boundary method, as we have implemented it, is formally second-order in space and first order in 
time, but, for problems with sharp interfaces that do not have smooth solutions, it is limited to first-order 
accuracy in space and time. Thus for our problem we expect only first order accuracy. The convergence 
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rates in time given in Table [T] on page |23| show first-order convergence in time as is expected. In Figure [TO] 
a), we depict the exact values of Ep{q{T); At) for g = X and q = u. We show log2 in the x 
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on page 

and y axes so that the empirical convergence rates from Table [T] on page |23] appear as the slope of the line 
segments. 

4-.3.3. Empirical Estimate of Convergence Rate in Space 
For this refinement study, we define a measure of error by 



E,iqinh)^\\q\T)-ll'^ (<Z"/'(T)) 



(39) 



which is the error difference at time t = T in a computed quantity, q, using a spatial refinement of a half. 
In this definition, if^^ is the restriction operator from a fine to a coarse grid. Then, an empirical estimate 
for the convergence rate is calculated using 



(40) 

We note that the estimates for convergence rates given by ( p38| and ( |40| ) have a fairly simple derivation using 
a Taylor series expansion (see [13] or |21jV 

In the spatial refinement analysis, we did not refine the Lagrangian grid with the Eulerian grid, so 
the same number of Lagrangian points were present in all of the simulations. In addition, full weighting 



restriction is used in the definition of the error, (39 1, for the error in u. We also used a fixed timestep of 
At = 10^^ until t = T = 0.01s for all of these simulations. The computed convergence rates from this 
refinement are provided in Table |2] on page |23| The oo-norm convergence rates given in Table [2] on page 
[23] show greater than first-order convergence in space for the error in the Lagrangian variable X and in the 
Eulerian variable u. The seemingly large convergence rates for the lower resolution grids {h — jq, -^21 'ii 
) can be explained by the fact that using a; = in the Dirac delta approximations does not allow the 
Lagrangian forces to be adequately represented in the Eulerian grid. This leads to larger errors in the 
coarse-grid simulations. Therefore, the best estimates for the convergence rates are the ones using the three 
resolutions all obeying uj > h given in the 4"^ and 7^^ columns of Table [2] on page 23 In Figure 10 on page 
b), we depict the exact values of Ep{q{T); h) for q = X and q = u. We show log2 in the x and y axes so 
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that the empirical convergence rates from Table |2] on page [23] appear as the slope of the line segments. We 
discuss possibilities for improvement in the convergence rates later in the conclusion sections. 

We also used a grid refinement analysis to find the empirical convergence rate with spatial refinement 
when the density of the biofilm is two times that of the surrounding fluid. This analysis was done to show 
that the first order convergence rate is maintained with the increased density in the biofilms. The results of 
this convergence analysis are shown in Table [3| on page 24 and Figure 10 on page 23 c). 

Finally, we compute the empirical convergence rates for our 2D simulation with variable viscosity. In 
this convergence study, (Table [4| on page k4^ and Figure [TO^ on page [23{^d) ) we use a non-dimensionalizec 
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value of biofilni viscosity of fimax = 500, which means that the viscosity at the location of a Lagrangian 
node is 500 times that of the surrounding fluid. First-order convergence in space is maintained, even with 
this very large biofilm viscosity. 




Figure 10: Empirical error estimates with (a) temporal refinement: Ep(q(T); At) is the p-norm of the error as defined by l |37| , 
(b) spatial refinement with constant density and viscosity: Ep{q{T); h) is the p-norm of the error as defined by | |39[ , (c) spatial 
refinement and increased biofilm density, (d) spatial refinement and increased biofilm viscosity. In all plots, we show log2 in 
the X and y axes so that the empirical convergence rate appears as the slope of the line segments. 
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Table 1: Empirical convergence rates with temporal refinement. rp{q(T); dt) is the convergence rate in the variable, q, a,t t = T 
using the p-norm and the three time steps dt, dt/4. 
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Table 2: Empirical convergence rates with spatial refinement. rp{q{T); h) is the convergence rate in the variable, q, at t = T 
using the p-norm and the three Eulerian step sizes h, '>/2, '»/4. 
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r2{q{T); 


^2(g(T);g^) 


^oo(g(r); j^) 


^oo(g(T); ^) 
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3.42 


1.19 
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1.20 
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1.83 
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1.88 



Table 3: Empirical convergence rates with spatial refinement and increased biofilm density. rp{q(T) \ h) is the convergence rate 
in the variable, q, at t = T using the p-norm and the three Eulerian step sizes h, ''/2, '>/4. In this experiment, the density of 
the biofilm is double that of the surrounding fluid. 
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r2{q{T);^) 
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roc{q{T): ^) 
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1.04 


0.78 


1.67 


0.86 
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2.04 


1.29 


1.82 


1.31 



Table 4: Empirical convergence rates with spatial refinement and increased biofilm viscosity. rp{q{T); h) is the convergence 
rate in the variable, q, a,t t = T using the p-norm and the three Eulerian step sizes h, '1/2, '1/4. In this experiment the viscosity 
of the biofilm is 500 times that of the surrounding fluid. 

4- .3. 4- Time-Step Stability Restrictions 

Finally, we investigated the stability of the method computationally as it depends on the spatial and 
temporal refinement and the stiffness of the springs. Analytically, stability applies to a numerical scheme 
and not to a computational run, but here we follow Mori and Peskin in [25 and give a simple definition of the 



stability for each computational run. Using the square of the 2-norm defined by (36l on u (i.e. ||u||j^) gives 
a value which is proportional to the kinetic energy in the system. We call the simulation stable if magnitude 
of the total velocity (as measured by the total kinetic energy) does not have a time of extreme growth during 
the simulation. Moreover, this kinetic energy should remain relatively close to the value of the total kinetic 
energy in the case of no biofilm. Using this definition of stability, we found, through experimentation with 
many combinations of h, At, and F^axi that we have timestep restrictions that scale with the mesh-width, 
h, and with the maximum Lagrangian force, Fmax- The restrictions are approximately given by 



At < Cih 



and 



C2 



F„ 



where Ci and C2 are positive proportionality constants. Specific values of Ci and C2 change depending on 
the parameters of the simulation. In future simulations, we hope to avoid these timestep restrictions by 
using an implicit or semi-implicit method as is done in [25 and [26 . All of the simulations shown in this 
work and used in the convergence testing used time-steps satisfying these two restrictions. 

4-.4- Three- Dimensional Validation 

In this subsection, we provide some numerical evidence validating the 3D simulations. We first validate 
in the absence of a biofilm using the exact laminar flow solution. Then, we validate the multigrid method 
in the presence of a biofilm and, finally, we provide the empirical convergence rates for the simulation in the 
presence of a biofilm. 
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We first tested the rate of convergence of our method on the laminar flow case without the interference 
of a biofilm. To illustrate the convergence rate in the absence of a biofilm, we started with an initial velocity 



profile that is one-half of that of the laminar flow velocity profile, given by (34 1. We ran the simulation 



enough timesteps until the approximate solution converged, with only discretization error remaining, to the 
exact solution for six spatial step sizes, h — {|, |, j^, g^, j^}- We computed the discretization error 



(using the exact laminar solution, (34 1, for computations) for each of the step sizes and found that the error 
is 0{h^). This can be seen in Figure 11 on page 25 where on the vertical axis we have the log2 of the error 



so that the convergence rate appears as the slope in the plot. 

Next, in Figure [9] on page 21 |b), we depict (for the pressure computation) the work units required to 
reach the minimum residual error as a function of allowed levels in the multigrid approach. This again 
implies that the multigrid method correctly accelerates the convergence of our iterative method for the 3D 
simulations with a biofilm. Note that there is only a slight reduction in the number of required work units 
with the addition of a 6"^ level in the multigrid, and we saw no reduction with 7 levels, so we use at most 6 
levels in our 3D solvers. 
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Figure 11: We show the convergence rate in the laminar flow case. These points are discretization errors for each of the spatial 
step sizes ( ^ , g , , i 64 ' T28 ) ' 



Finally, as was done for the 2D case in |4.3.3| we compute the empirical convergence rates for our 3D 
simulation in the presence of the biofilm shown in Figure[7]on page 19 h) with all of the same fiuid parameters 



used in the 2D analysis. Using the p-norms defined above, we can compute the convergence rates using (38 1 



and ( [40| (see Figure 12 on page 26 'a) and Table [s] on page 26). For the temporal convergence analysis, we 
used uj = and h — This analysis resulted in first-order convergence in all measures except rp{q{T); At) 
in which it has an average convergence rate of about 0.6. Next, we found empirical convergence rates for 
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spatial refinement (see Table [6] on page 27 and Figure 12 on page 26 b)). As expected, we observe a greater 
than first-order convergence rate in both the Eulerian velocity, u, and the Lagrangian position, X. Next, 
we conducted a spatial refinement analysis with a biofilm that has double the density of the surrounding 
fluid (see Table [t] on page |27| and Figure [T2| on page |26[ c)). Finally, we did the spatial refinement study for 
our simulations with firnax = 500, and again achieved first-order spatial convergence (see Table [s] on page 



27 and Figure 12 on page 26 d)) 
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Figure 12: Empirical errors in the 3D simulations with (a) temporal refinement: Ep{q{T); At) is the p-norm of the error as 
defined by | |37[ l, (b) spatial refinement with constant density and viscosity: Ep(q{T); h) is the p-norm of the error as defined 
by | |39[ , (c) spatial refinement and increased biofilm density, (d) spatial refinement and increased biofilm viscosity. We show 
log2 in the x and y axes so that the empirical convergence rate appears as the slope of the line segments. 
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r2(g(T);.0004) 


r2(g(T);.0002) 


r2(g(r);10-4) 


Too (g(r);. 0004) 


r^(<z(T);.0002) 


roo(g(T);lG-4) 


u 


0.71 


0.35 


0.64 


1.04 


1.18 


0.83 


X 


0.90 


1.03 


1.11 


0.90 


1.02 


0.88 



Table 5: Empirical convergence rates in the 3D simulations with temporal refinement are shown for u and X. rp(g(T); At) is 
the convergence rate in the variable, q, at t = T using the p-norm and the three Eulerian step sizes At, ^*/2, 
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q 




r2{q{T); 


^2(g(T);g^) 


^oo(g(r); j^) 


^oo(g(T); ^) 




u 


2.09 


1.44 


1.02 


2.89 


2.11 


1.47 


X 


3.12 


1.57 


1.47 


2.70 


1.33 


1.23 



Table 6: Empirical convergence rates in the 3D simulations with spatial refinement are shown for u and X. rp(q{T); h) is the 
convergence rate in the variable, q, at t = T using the p-norm and the three Eulerian step sizes h, h/2, h/i. 



q 




r2{q{T); ^) 


r2{q{Th^,) 


roo{q{T); j^) 


roo{q{T); ^) 


roo{q{T)\ ^) 
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1.95 


1.52 


1.40 


2.52 


1.61 


1.73 


X 


1.7 


1.28 


1.55 


1.29 


1.40 


1.56 



Table 7: Empirical convergence rates with 3D spatial refinement and increased biofilm density. rp(g(T); h) is the convergence 
rate in the variable, g, at t = T using the p-norm and the three Eulerian step sizes h, ''/2, In this experiment, the density 
of the biofilm is double that of the surrounding fiuid. 
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r2{q{TYl,) 
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0.86 


0.57 


1.04 


0.92 


X 


1.47 


1.07 


1.48 


1.09 



Table 8: Empirical convergence rates with 3D spatial refinement and increased biofilm viscosity. rp{q(T); h) is the convergence 
rate in the variable, g, at t = T using the p-norm and the three Eulerian step sizes h, '»/2, h/i. In this experiment, the viscosity 
of the biofilm is 500 times that of the surrounding fluid. 



This concludes our validation section, and we now present the results of our numerical simulations. 

5. Simulations Results 

In this section, we present the results of our numerical simulations. First, we briefly discuss the reality 



of elastic forces in biofilms. Then we provide 2D results in 5.2 and 3D results in 5.3 

5.1. Discussion of Elastic Maximum Force, F^ax 

We now provide a brief discussion of the physical reality of the values of F^ax used in our simulations. 
The cohesive strengtl^Ti biofilms has been found experimentally to be highly heterogeneous, with repeated 
experimental measurements on the same biofilm yielding vastly different strength measurements. For ex- 
ample, 49 cohesive strength measurements taken on only two samples of Staphylococcus epidermidis yielded 
measurements between 61-5182 Pa fT. These biofilms were grown on a 22 mm diameter disc rotating at 
75 ^°t/min so the fastest speed, ^ 86m™/s, was at the perimeter of the disc (i.e. very slow flow growth 
conditions). The ad/iesiu^^nd cohesive strengths have also been shown to vary significantly with changes 
in growth conditions such as flow rate and nutrient concentration. Changes in these growth conditions 
influence the amount of ECM production in the biofilm as well as the compactness of the biofilm, which 
has a direct effect on its strengths |10l [271 [3] . We note here that the required values we find for F„iax for 



^The cohesive strength is a measure of the forces that interconnect the biofilm's cells. 
®The adhesive strength is a measure of the forces that connect a biofilm to the surface. 
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Parameter Values for the Simulations 


Tube Radius 


25 X 10"^ m 


Fluid Dynamic Viscosity 


1.0 X 10-'^kg/,„., 


Fluid Density 


998 kg/m-' 


Maximum Fluid Velocity 


10"^ m/s 



Table 9: The values of parameters used in the 2D simulations. 

the biofilms to remain attached in our 2D and 3D simulations are consistent with the cohesive strength 
measurements provided in jT|. Since the diameter Staphylococcus epidermidis is about 1 /im. in 3D. we 
multiply the cohesive strengths by 1 /im^ to get an approximation for the range of forces on the surface area 
of one cell. Using the range of 61-5182 Pa yields a range of forces from 6.1 x 10^^^ N to 5.18 x 10^^ N. 
In 2D, we multiply the cohesive strengths by the cell diameter to get a rough approximation for the range 
of forces on the surface perimeter surrounding one cell. Using the range of 61-5182 Pa yields a range of 
forces from 6.1 x 10~^ N to 5.18 x 10~^ N. Our values for Fmax are at the low end of these ranges. The 
actual strength of the biofilm is most likely larger than our F„iax values since the positional data was from 
a biofilm that was not fragmenting in the flow conditions in which it was grown, and we used the same flow 
conditions our simulations. Thus, in order to see detachment under these flow conditions, we had to lower 
the value of Fmax- We could alternatively increase the flow rate to necessitate a larger Fmax requirement 
to avoid detachment. One eventual goal of this work is that, if the approximate value of F^ax is known for 
a particular type of biofilm, then our simulations can be used to predict the fiow rates required to break 
different shaped biofilms. 

5.2. Two-Dimensional Simulations 

In this section, we provide results from our 2D simulations, which represent a cross-section of a biofilm 
attached to the inside of a tube and subjected to fiuid fiow in a computational domain of 150 /xm by 50 /xm. 
The parameters for our simulations are given in Table |9] on page [28] 

In all simulations, we implement a breaking condition on the springs of two times the rest length. The 
initial configuration for the biofilm in these simulations is shown in Figure [7] on page |19[ The spring 
connections between Lagrangian nodes are put in place at the beginning of the simulation with every node 
connected to every other node less than dc away (the reason for this connection distance is given above in 



4.2 1. The mushroom shaped biofilm has a height of about 8.5 /im and width of about 8 /im (width of about 
2 /im at the thinnest part) . We use a non-dimensionalized lo = to match the radius of Staphylococcus 
epidermidis and choose h — in all of the simulations shown, so that uj > h. 

In the first simulation, the maximum spring force, Fmax, is set to 5.00 x lO^^N, and the results are 
provided in Figure [13] on page |30[ The biofilm bends over in the fiow, and the connections in the thin part 



of the biofilm break as they stretch too far. The streamlines in (b), (c), and (d) of Figure 13 on page 30 and 
in all of the other 2D simulation plots follow the trajectories given by the velocity field, u. 

28 



We point out that the values of the spring constants are well within physically realistic values (see the 



discussion in 5.11, although, in this work, we have chosen these values for the qualities they give to the 
simulations rather than experimental evidence of the elastic strength of biofilms. For example, in these 2D 
simulations, we investigated several simulation runs with various spring constants until we obtained those 
that exhibited the above described behaviors. 

Next, we conducted a simulation of the same mushroom shaped biofilm with all of the same parameter 
values, but we gave the biofilm additional density of — 998kg/m^ compared to the ambient fiuid. We 
know this density is larger than what is seen in actual biofilms (at most 20% greater density than water 
|241 137|). but we chose it to show an exaggerated example of increasing the biofilm density. These result is 
provided in Figure 14 on page 31 b) and illustrates that the added density essentially adds momentum to 



the biofilm. This additional momentum causes the biofilm to curl over into the slower flow region and thus 
prevents detachment. 

Finally, we conducted a simulation of the same mushroom shaped biofilm with all of the same parameter 
values as the first simulation, but we increased F^ax to 5.00 x 10^^ N. The effect of these stronger springs 
is that the thin part of the biofilm does not stretch enough to break the connections. The result is depicted 



in Figure 14 on page 31 'c). We can see from these simulations that either increasing the biofilm density or 
strengthening the springs causes similar results, but, with the increased density, the biofilm has more of a 
curling action. 

In the next simulation, we use all of the same parameters described in the first simulation, with the 
addition that the biofilm has a 500x larger viscosity than the surrounding fiuid, so [imax = 0.5kg/ms. 
Comparing simulation results illustrated in Figure 



14 



on page 



31 ^d) and Figure 



14 



on page 



31a), which 



show the biofilm configurations just before detachment, we observe a longer time until detachment in the 
high viscosity case. This is the expected outcome of increasing the viscosity in the biofilm. We note here 



that we used lo = in the equation for /i(x) in (22) because we wanted to spread additional viscosity over 
the same region that the elastic forces are spread to. We achieve an even longer detachment time in the 



simulation by widening the infiuence of additional viscosity by using, for example cj 



50- 
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Figure 13: 2D Simulation of a mushroom shaped biofilm with the same density as the surrounding fluid. Time is in seconds 
and the distance is in microns. In this simulation, po = 998kg/m^, p;, = 0, Fmax = 5.00 X 10~^N. The streamlines follow 
the velocity field. In this simulation, the top of the biofilm stretches in the flow, and the top breaks off as the connections in 
the the middle separate as they exceed the breaking criteria of twice the rest length. As expected in a laminar shear flow, the 
broken piece then tumbles end over end through the flow. 
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(a) 



(b) 



(c) 
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Figure 14: Snapshots of cell distributions resulting from 2D simulation of mushroom shaped biofilm with four different property 
configurations. In (a), the biofilm has the same density and viscosity as the surrounding fluid (just before detachment), (b) 
increased Fmax, (c) density twice as large as the surrounding fluid density, and (d) viscosity is 500 X that of the surrounding 
fluid (just before detachment). The streamlines follow the velocity field. 



5.3. Three- Dimensional Simulations 

In this section, we provide results from our 3D simulations, which use the same parameter values as 
in the two dimensional simulations (see Table |9] on page 28 ) . The difference is that the simulation in 3D 
represents flow through a square shaped tube with a side length of 50 /im. Note that these 3D simulations 
reproduce qualitatively the same results as in the 2D ones. 

Our 3D simulations were run on a 50 x 50 x 150 /im computational domain. We simulate on a mushroom 
shaped biofilm with a height of about 8.5 /im and a diameter of about 7.5 /im. This shape is carved from 
the same set of data points described in |4.2[ The spring connections between Lagrangian nodes are put in 
place at the beginning of the simulation, with every node connected to every other node less than dc = 3 /im 
away. Note that, for the 3D simulations, we increased dc slightly to establish enough connections in the 
biofilm. We again use u = to match the radius of Staphylococcus epidermidis and choose h = in 
all of the simulations shown, so that uj > h. In the first simulation, the maximum spring force, Fmax, is 
set to 1.25 X 10^^^ N. We again chose the value of these spring constants in order to illustrate specific 
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behaviors. The results of the first simulation are shown in Figure [15] on page |33[ The mushroom shaped 
biofilm bends over and stretches in the flow. The connections in the midsection of the biofilm exceed their 
breaking length and the top of the biofilm breaks off into the flow. Next, we ran a simulation of the same 
mushroom shaped biofilm, but we added ph = 998 kg/m^ additional density to the biofilm compared to the 



surrounding fiuid. The final result is provided in Figure 16 on page 34 'b) and illustrates that the added 
density increases the momentum of the biofilm. This allows for the mushroom to curl over into the fiow and 
increases the time until detachment. We also ran a simulation of the same mushroom shaped biofilm, but we 
increased F^ax to 1 x 10^^^ N and kept the biofilm density the same as the surrounding fiuid. The result is 



provided in Figure 16 on page |34[ c). The effect of these stronger springs is that the thin part of the biofilm 
does not stretch enough to break the connections. We can see from these simulations that either increasing 
the biofilm density or strengthening the springs causes similar results, but with the increased density the 
biofilm just curls over. 

Finally, we provide one 3D simulation to show that, with increased biofilm viscosity, they produce 
qualitatively the same behavior as in the 2D case. In the simulation result shown in Figure [16] on page 



34 d), we use the same parameters as the first 3D simulation, but we use a viscosity in the biofilm that is a 



factor of 500 times that of the surrounding fiuid. Just as in the 2D case, this results in a longer time until 

on page 



detachment (compare time in Figure 16 on page 34 ^a) and Figure 16 



34;d)). 



For ease of comparison, we now provide the figures in the order in which they were discussed in this 
section. 
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(c) X z (d) ^ 

Figure 15: Full 3D simulation of a mushroom shaped biofilm with the same density as the surrounding fluid. The time is in 
seconds and the distance is in microns. As the biofilm stretches in the flow, the strain in the midsection exceeds the breaking 
length of the connections, and the top of the biofilm breaks off into the flow. Then the broken piece tumbles end over end 
through the flow, and the base retracts back. In this simulation, po = 998kg/m^, pj, = 0, F^ax = 1-25 X 10~^^ N. 
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t=0.12 



(a) 




(b) 






Figure 16: Snapshots of cell distributions resulting from 3D simulation of mushroom shaped biofilm with four different property 
configurations. In (a), the biofilm has the same density and viscosity as the surrounding fluid (just before detachment), (b) 
increased Fmax, (c) density twice as large as the surrounding fluid density, and (d) viscosity is 500 X that of the surrounding 
fluid (just before detachment). 



6. Realistically Shaped Biofilm Simulation 



The biofilm shapes used in |5 . 2| and |5 . 3| were intentionally carved from the data in a way to provide a weak 
point at which it would be most likely to break. This was done in order to illustrate the effects of varying 
the different parameters in the simulation. In this section, we provide results of the simulation on a biofilm 
that is a subset of points taken directly from the real biofilm data set. In reality, this biofilm was surrounded 
by more cells on all sides, which would change the behavior of the fluid structure interactions. However, we 
use this to show the results of the simulation on a real top heavy biofilm shape that was grown in a lab. 
The Staphylococcus epidermidis data set discussed in |4.2| was supplied as positions in three 30 x 30 x 15 /im 
sub-domains of a biofilm. 



In Figure 17 



on page 



36 a), we show the biofilm taken from a 2 x 30 x 15/im subset of one data set that 



has been connected with dc = 2.8 /im. For the 2D representation, we collapse the 2 /im dimension, leaving 
only the (x, y) coordinates of the data. The most interesting feature of this biofilm is that in the region 
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from X = 60 /im to x = 67 /im the biofilni exhibits a mushroom shape similar to the one we used in |5.2[ 
We now provide two simulation results on this realistically shaped biofilm. The first simulation (see 
36) uses a biofilm density equal to the surrounding fluid and uses F^ax — 7.5 x 10^^ N. 



Figure 



17 



on page 



In this simulation, the mushroom shaped part pushes against the biofilm behind it. then rolls over the top of 
it as it breaks from its base, forming a long streamer-like biofilmQlhen the streamer breaks completely off 
leaving two distinct attached structures. In the second simulation, we use a biofilm density of pf, = 120 '^s/m^ 
and kept everything else the same. Although the density is only 12% larger than the surrounding fiuid, it 
has a large impact on the outcome of the simulation. In this simulation, the effect of the increased density 

on page 



of the biofilm is a longer breaking time (compare time in Figure 18 



37 'b) and Figure 18 



on page 



p7|c)). This occurs since the increased momentum causes the first detached piece to continue further down, 
pulling the whole streamer lower (compare the height of the detaching pieces). The fluid forces continue to 
push the streamer until it breaks into the flow. 

In the final simulation, we show that increasing the viscosity in the "realistically" shaped biofilm has 
larger impact on the results than in the case of the previous standalone mushroom shaped biofilm. In 
this simulation, we increased the biofilm viscosity to 50 x the surrounding fiuid with fimax — 0.05kg/m s. 
Although this is 10 x less than in the previous variable viscosity simulation, it has a larger impact on this 
wider biofilm, doubling the detachment time from the case of constant viscosity (compare Figure 18 on page 
37 b) and Figure 18 on page [37|d)). 



Streamers are a natural occurrence in biofilms. Examples in 1171 1391 . 
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Figure 17: Simulation on a 2D slice of a real biofilm with the same density as the surrounding fluid. Time is in seconds and the 
distance is in microns. In this simulation, po = 998kg/m^, pi, = 0, Fmax = 7.5 X 10~^ N. The streamlines follow the velocity 
field. In this simulation, the mushroom shaped part pushes against the biofilm behind it (b), then rolls over the top of it as it 
breaks from its base (d). Then a large portion of the biofilm breaks completely off leaving 2 distinct bases (f). 
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Figure 18: Snapshots in time showing the detachment time and configuration at detachment of a 2D slice of a real biofilm with 
initial configuration in (a). In (b), biofilm has the same density and viscosity as the surrounding fluid, (c) 12% larger density, 
and (d) 50 X larger viscosity than the surrounding fluid. The streamlines follow the velocity field. 

6.1. Numerical Concerns 

We note that there remains one pressing concern that will be the focus of our future work. When adapting 
the multigrid scheme for large values of ^max, we do not achieve expected speed ups in convergence rates. 
We use restriction to transfer the viscosity to the coarse grids works for small values of iJimax , but we found 
that, for larger values of ^max, this technique leads to a very slowly converging solver. Intriguingly, we found 
that using our restriction operator to define the coarse grid viscous values and then scaling the values leads 
to faster convergence. Specifically, we define the coarse grid viscosity as 

= lihlf^f^^hi^), (41) 
where Ih denotes the grid whose mesh width is I times h {I — 2, 4, 8, 16, . . .), jih G (0, 1] is the scaling which 



maximizes the convergence of the solver, and /i/i(x) is defined by (22 1. Through repeated experimentation 



we found that using '^ih G [.7, 1] resulted in the fastest convergence rates in the 2D and 3D simulations. 
This approach is admittedly ad-hoc. However, the consistency with which we achieved dramatic speed ups 
strongly suggests the existence of an underlying mathematical principle to be discovered. 
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For our future work, the highest priority is to resolve the problem of slow convergence for large ^max- 
There are three approaches that may lead to resolving this issue. The first would be to mathematically derive 
optimal values for the scaling parameters, 7;^. Second, we could ensure that our discretization satisfies the 
Galerkin condition. Lastly, and probably the best choice, would be to re-implement the geometric multigrid 
as an algebraic multigrid method (AMG, see Ch. 8 of [7]). 

7. Conclusions 

In this work we developed a simulation to model the flow-induced fragmentation of biofilms. In this 
simulation, we have provided a way to adjust the biofilm density and viscosity, which had not been addressed 
in previous IBM biofilm models. We also have control of the fluid flow rate, density, viscosity, and elastic 
forces within the biofilm. We used experimentally measured biofilm bacterial cell locations as initial positions 
for our Lagrangian nodes. This is dramatically different than the traditional IBM, in which methods usually 
refine the Lagrangian mesh along with the Eulerian mesh. We adapted the Dirac delta approximation to 
scale with the radius of the bacteria rather than with the mesh width. This implies that the information that 
transfers from the Lagrangian grid to the Eulerian grid (i.e., density, viscosity, and elastic force) is spread 
over a set distance rather than scaling by the mesh width, h. This adapted Dirac Delta approximation 
improves our numerical convergence rates as well. 

We used a projection method to split the incompressible Navier-Stokes equations to solve separately for 
an intermediate velocity and the pressure, and then used a Gauss-Seidel iterative method with multigrid 
to solve the resulting equations. Using an iterative solver, as opposed to a spectral method, to solve these 
systems was necessitated by the fact that biofilms have spatially varying density and viscosity. With this 
solver we achieved first order convergence in both space and time. 

For the numerical simulations, we carved a mushroom shaped biofilm from the bacterial cell locations 
and ran simulations with varying parameters. We first ran the simulation on a simplistic shape in order 
to validate the effect of the various parameter changes on the outcome of the biofilm. By adjusting the 
maximum elastic force, F^ax in the biofilm, we controlled the detachment phenomenon. We also showed 
that slight changes in the density of the biofilm has a large effect on the outcome of the simulation. This is an 
important conclusion as usually modelers ignore the differences in biofilm density. Finally, we showed that 
we can increase the detachment time in the simulations by increasing the viscosity of the biofilm. Finally, 
we ran simulations on more realistically shaped biofilms, which showed how a larger biofilm with different 
shapes will react to fluid flow forces. Adjusting these parameters will be a necessary component when we 
attempt to match these simulations to experimental data. 
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8. Future Work 



There are several directions in which we plan to take this research in the future. For example, there 
are straightforward ways to include more biologically realistic terms to interpret biofilm internal stress 
dynamics, cell volume, and fragmentation dynamics. Additionally, there are approaches that may greatly 
improve our numerical scheme convergence and stability, including alternative multigrid algorithms, implicit 
discretizations, and improved immersed boundary implementations. 

In its current form, our simulations could be used to make predictions in detachment times of biofilms, 
as well as general behavioral responses of biofilms to various flow conditions. We plan to work closely 
with experimentalists to formulate accurate viscoelastic models for the biofilms and modify our constitutive 
equations for stress and elasticity to account for these model choices. We also plan to include the fact that 
bacterial cells displace fluid. While we have adapted the Dirac delta function approximation to transfer 
the cell parameters (F, p, fi) to the Eulerian grid, the current simulation does not actually assign a size to 
the cells. As a result, the cells are free to pass through each other. We flrst plan to alter the model for 
the bacterial cell so that it displaces fluid. An important step in this process will be to identify a collision 
detection strategy. We will base ours on potentials for electrostatic, steric, and Van der Waals forces. 

Our current simulation uses a spring-breaking criteria of double the rest length, and we also assume that 
the bonds are linearly elastic until the breaking point. This is not an accurate assumption, as it is known 
that biofllms are composed of polymer based ECMs. These structures are linearly elastic for small strains 
and then experience plastic deformation (permanently altering the bonds in the ECM and thus the rest 
length) before finally fracturing. In the future, we will use biofilm yielding data from experiments such as 
[IJ to determine accurate approximations for yield points and fracture points in the biofilm. We plan to 



include plasticity into the simulations by changing the equations for stress, (20 1, when the bond has been 
stretched beyond its yield point. 

We have several plans for improving our numerical method. Our current simulation is limited to first- 
order accuracy. Guided by the results in [8\, we will accurately derive the numerical boundary conditions 
for our projection method to ensure second order accuracy for both velocity and pressure computations. 
To improve the accuracy of the immersed boundary method, we could also adapt our modeling method 
to either an immersed interface method ([22J) in which we adapt the finite difference approximations close 
to the interface or a blob projection immersed boundary method as discussed in [11 in order to obtain 
second-order spatial accuracy. Another limitation of our current numerical scheme is the time-step stability 
restrictions, which limit the size of the elastic forces between the cells. We plan to eliminate these restrictions 
altogether by changing to a semi-implicit or implicit method of transferring the data between the Eulerian 
and Lagrangian grids, as is shown by Newren, Fogelson et al., in [26J. 

Finally, with large biofilm densities and viscosities, our multigrid method in its current formulation does 
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not converge as fast as expected. We plan to fix this by appropriately adapting our implementation of the 
geometric multigrid (by satisfying the Galerkin condition) or by changing to an algebraic multigrid approach. 
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Appendix A. 

In this appendix we provide a list of variables and parameters used in this paper. 
b Dashpot Damping Coefficient 
do Average Spring Rest Length 
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5 Dirac Delta Function 

Discretized Dirac Delta Function from Peskin 

6 Our Modified Discretized Dirac Delta Function 

D Spatial dimension, D = 2 for 2D simulations and = 3 for 3D simulations 

ej Unit Vector in the i*'^ direction 

7] Total Number of Lagrangian Points 

f Eulerian Force Density 

F Lagrangian Force 

Fmax Maximum Lagrangian Force 

h Spatial Discretization of finest grid 

K Hookean Spring Coefficient 

/U Dynamic Fluid Viscosity 

IJ-max Maximum Biofilm Viscosity 

p Pressure 

q = {q, r, s) Lagrangian Coordinates 
p Density 

Pq Uniform Fluid Density 

Pb Additional Density in Biofilm 

s Lagrangian Node Marker 
t Time 

T Tension in Spring 

u Eulerian Velocity 

u Intermediate Velocity 

U Lagrangian Velocity 

X = (xi,a;2,a;3) Cartesian Coordinates 
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Appendix B. Scaling Parameters 



Scaling Parameters 



Description 


ScaUng 
parameter 


Primary 
dimensions 


Specific values 
chosen for 

simulations 


Characteristic Length 


L 


{L} 


50 microns 


Characteristic Speed 


Uq 




10-^™/s 


Characteristic Frequency 


T 


{t} 


Is 


Reference Pressure Difference 


PO - PLt^be 




.8144 Pa 


Characteristic Density 


Po 




998 kg/m^ 


Characteristic Viscosity 


/" 


{mL-H-'} 




Characteristic Force Density 


/o 




varies 



Table B.IO: Shows the scaling parameters and their descriptions. 
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